/*
 * TurtleKit - An Artificial Life Simulation Platform
 * Copyright (C) 2000-2010 Fabien Michel, Gregory Beurier
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 */
package edu.turtlekit2.tools.fillers;

import java.util.Random;


/**
 * A Perlin noise generator to generate noise of grid variables.
 * @author G. Beurier
 * @version 1.0 - 6/2008
 * Adapted from source code by Carl Burke. Adapted from copyrighted source code by Ken Perlin and F. Kenton Musgrave to accompany: Texturing and Modeling: A Procedural Approach Ebert, D., Musgrave, K., Peachey, P., Perlin, K., and Worley, S. AP Professional, September, 1994. ISBN 0-12-228760-6
 */
public class PerlinSolidNoiseGenerator implements SolidNoiseGenerator
{
    // *** METHODS OF TERRAIN GENERATION CURRENTLY SUPPORTED

    static public final int METHOD_BASIC = 1;
    static public final int METHOD_MULTIFRACTAL = 2;
    static public final int METHOD_HETERO_TERRAIN = 3;
    static public final int METHOD_HYBRID_MULTIFRACTAL = 4;
    static public final int METHOD_RIDGED_MULTIFRACTAL = 5;

    // *** PRIVATE DATA TO DRIVE TERRAIN CALCULATIONS

    private int method;
    private double H;
    private double lacunarity;
    private double octaves;
    private double offset;
    private double gain;
    private double[] point;

    private Random rgen;

    public PerlinSolidNoiseGenerator()
    {
        rgen = new Random();
	init_noise();
	point = new double[3];
	method = METHOD_BASIC;
	H = 0.5;
	lacunarity = 2.0;
	octaves = 7.0;
    }

    public PerlinSolidNoiseGenerator(double hIn, double lacIn, double octIn)
    {
        rgen = new Random();
	init_noise();
	point = new double[3];
	method = METHOD_BASIC;
	H = hIn;
	lacunarity = lacIn;
	octaves = octIn;
    }

    public PerlinSolidNoiseGenerator(int methIn, double hIn, double lacIn,
	double octIn, double offIn)
    {
        rgen = new Random();
	init_noise();
	point = new double[3];
	switch (methIn)
	{
	  case METHOD_MULTIFRACTAL:
	    method = METHOD_MULTIFRACTAL;
	    H = hIn;
	    lacunarity = lacIn;
	    octaves = octIn;
	    offset = offIn;
	    break;
	  case METHOD_HETERO_TERRAIN:
	    method = METHOD_HETERO_TERRAIN;
	    H = hIn;
	    lacunarity = lacIn;
	    octaves = octIn;
	    offset = offIn;
	    break;
	  case METHOD_HYBRID_MULTIFRACTAL:
	    method = METHOD_HYBRID_MULTIFRACTAL;
	    H = hIn;
	    lacunarity = lacIn;
	    octaves = octIn;
	    offset = offIn;
	    break;
	  default:	// don't know which method, so do basic
	    method = METHOD_BASIC;
	    H = hIn;
	    lacunarity = lacIn;
	    octaves = octIn;
	}
    }

    public PerlinSolidNoiseGenerator(double hIn, double lacIn,
	double octIn, double offIn, double gainIn)
    {
        rgen = new Random();
	init_noise();
	point = new double[3];
	method = METHOD_RIDGED_MULTIFRACTAL;
	H = hIn;
	lacunarity = lacIn;
	octaves = octIn;
	offset = offIn;
	gain = gainIn;
    }

    ///** RANDOM NUMBER INITIALIZATION AND GENERATION ***

    private double rseed;

    public void setSeed(double s)
    {
	rseed = s;
	rgen.setSeed(Double.doubleToLongBits(rseed));
	init_noise();
    }
    public double getSeed() {return rseed;}

    /************************************************************
     * Methods specific to noise synthetic terrain generation
     ************************************************************/
    /************************************************************
     * Supporting/filtering methods
     ************************************************************/

    public double bias(double a, double b)
    {
	return Math.pow(a, Math.log(b) / Math.log(0.5));
    }
    
    public double gain(double a, double b)
    {
	double p = Math.log(1. - b) / Math.log(0.5);

	if (a < .001)
	    return 0.;
	else if (a > .999)
	    return 1.;
	if (a < 0.5)
	    return Math.pow(2 * a, p) / 2;
	else
	    return 1. - Math.pow(2 * (1. - a), p) / 2;
    }

    private double vec[];
    public double turbulence(double v[], double freq)
    {
	double t;

	if (vec==null) vec = new double[3];
	for (t = 0. ; freq >= 1. ; freq /= 2)
	{
	    vec[0] = freq * v[0];
	    vec[1] = freq * v[1];
	    vec[2] = freq * v[2];
	    t += Math.abs(noise3(vec)) / freq;
	}
	return t;
    }
    
    /************************************************************
     * Initialization
     ************************************************************/

    private void normalize2(double v[]) // v.length == 2
    {
	    double s;
    
	    s = Math.sqrt(v[0] * v[0] + v[1] * v[1]);
	    v[0] = v[0] / s;
	    v[1] = v[1] / s;
    }
    
    private void normalize3(double v[]) // v.length == 3
    {
	    double s;
    
	    s = Math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
	    v[0] = v[0] / s;
	    v[1] = v[1] / s;
	    v[2] = v[2] / s;
    }
    
    private void init_noise()
    {
	int i, j, k;

	p = new int[B + B + 2];
	g3 = new double[B + B + 2][3];
	g2 = new double[B + B + 2][2];
	g1 = new double[B + B + 2];

	for (i = 0 ; i < B ; i++)
	{
	    p[i] = i;

	    g1[i] = rgen.nextDouble() * 2.0 - 1.0;  // -1.0 to 1.0

	    for (j = 0 ; j < 2 ; j++)
		g2[i][j] = rgen.nextDouble() * 2.0 - 1.0;  // -1.0 to 1.0
	    normalize2(g2[i]);

	    for (j = 0 ; j < 3 ; j++)
		g3[i][j] = rgen.nextDouble() * 2.0 - 1.0;  // -1.0 to 1.0
	    normalize3(g3[i]);
	}

	while ((--i) > 0)
	{
	    j = (int)(rgen.nextDouble() * B);
	    k = p[i];
	    p[i] = p[j];
	    p[j] = k;
	}

	for (i = 0 ; i < B + 2 ; i++)
	{
	    p[B + i] = p[i];
	    g1[B + i] = g1[i];
	    for (j = 0 ; j < 2 ; j++)
		g2[B + i][j] = g2[i][j];
	    for (j = 0 ; j < 3 ; j++)
		g3[B + i][j] = g3[i][j];
	}
    }

    /************************************************************
     * Noise generation (interpolation) over 1,2, and 3 dimensions
     ************************************************************/

    /* noise functions over 1, 2, and 3 dimensions */
    
    static final int B = 0x100;
    static final int BM = 0xff;
    
    static final int N = 0x1000;
    static final int NP = 12;   /* 2^N */
    static final int NM = 0xfff;
    
    private int p[];
    private double g3[][];
    private double g2[][];
    private double g1[];
    
    public double s_curve(double t)
    {
	return t * t * (3. - 2. * t);
    }
    
    public double lerp(double t, double a, double b)
    {
	return a + t * (b - a);
    }
    
    /*
    #define setup(i,b0,b1,r0,r1)\
	t = vec[i] + N;\
	b0 = ((int)t) & BM;\
	b1 = (b0+1) & BM;\
	r0 = t - (int)t;\
	r1 = r0 - 1.;\
    */

    public double noise1(double arg)
    {
	int bx0, bx1;
	double rx0, rx1, sx, t, u, v;

	/* setup(0, bx0,bx1, rx0,rx1) */
	t = arg + N;
	bx0 = ((int)t) & BM;
	bx1 = (bx0+1) & BM;
	rx0 = t - (int)t;
	rx1 = rx0 - 1.;
	/***/

	sx = s_curve(rx0);

	u = rx0 * g1[ p[ bx0 ] ];
	v = rx1 * g1[ p[ bx1 ] ];

	return lerp(sx, u, v);
    }
    
    public double noise2(double vec[])  // vec.length == 2
    {
	int bx0, bx1, by0, by1, b00, b10, b01, b11;
	double rx0, rx1, ry0, ry1, q[], sx, sy, a, b, t, u, v;
	int i, j;

	/* setup(0, bx0,bx1, rx0,rx1) */
	t = vec[0] + N;
	bx0 = ((int)t) & BM;
	bx1 = (bx0+1) & BM;
	rx0 = t - (int)t;
	rx1 = rx0 - 1.;
	/***/
	/* setup(1, by0,by1, ry0,ry1) */
	t = vec[1] + N;
	by0 = ((int)t) & BM;
	by1 = (by0+1) & BM;
	ry0 = t - (int)t;
	ry1 = ry0 - 1.;
	/***/

	i = p[ bx0 ];
	j = p[ bx1 ];

	b00 = p[ i + by0 ];
	b10 = p[ j + by0 ];
	b01 = p[ i + by1 ];
	b11 = p[ j + by1 ];

	sx = s_curve(rx0);
	sy = s_curve(ry0);
    
	q = g2[ b00 ] ; u = ( rx0 * q[0] + ry0 * q[1] );
	q = g2[ b10 ] ; v = ( rx1 * q[0] + ry0 * q[1] );
	a = lerp(sx, u, v);

	q = g2[ b01 ] ; u = ( rx0 * q[0] + ry1 * q[1] );
	q = g2[ b11 ] ; v = ( rx1 * q[0] + ry1 * q[1] );
	b = lerp(sx, u, v);

	return lerp(sy, a, b);
    }
    
    public double noise3(double vec[])  // vec.length == 3
    {
	int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
	double rx0, rx1, ry0, ry1, rz0, rz1, q[], sy, sz, a, b, c, d, t, u, v;
	int i, j;

	/* setup(0, bx0,bx1, rx0,rx1) */
	t = vec[0] + N;
	bx0 = ((int)t) & BM;
	bx1 = (bx0+1) & BM;
	rx0 = t - (int)t;
	rx1 = rx0 - 1.;
	/***/
	/* setup(1, by0,by1, ry0,ry1) */
	t = vec[1] + N;
	by0 = ((int)t) & BM;
	by1 = (by0+1) & BM;
	ry0 = t - (int)t;
	ry1 = ry0 - 1.;
	/***/
	/* setup(2, bz0,bz1, rz0,rz1) */
	t = vec[2] + N;
	bz0 = ((int)t) & BM;
	bz1 = (bz0+1) & BM;
	rz0 = t - (int)t;
	rz1 = rz0 - 1.;
	/***/

	i = p[ bx0 ];
	j = p[ bx1 ];

	b00 = p[ i + by0 ];
	b10 = p[ j + by0 ];
	b01 = p[ i + by1 ];
	b11 = p[ j + by1 ];

	t  = s_curve(rx0);
	sy = s_curve(ry0);
	sz = s_curve(rz0);
    
	q = g3[ b00 + bz0 ] ; u = ( rx0 * q[0] + ry0 * q[1] + rz0 * q[2] );
	q = g3[ b10 + bz0 ] ; v = ( rx1 * q[0] + ry0 * q[1] + rz0 * q[2] );
	a = lerp(t, u, v);

	q = g3[ b01 + bz0 ] ; u = ( rx0 * q[0] + ry1 * q[1] + rz0 * q[2] );
	q = g3[ b11 + bz0 ] ; v = ( rx1 * q[0] + ry1 * q[1] + rz0 * q[2] );
	b = lerp(t, u, v);

	c = lerp(sy, a, b);

	q = g3[ b00 + bz1 ] ; u = ( rx0 * q[0] + ry0 * q[1] + rz1 * q[2] );
	q = g3[ b10 + bz1 ] ; v = ( rx1 * q[0] + ry0 * q[1] + rz1 * q[2] );
	a = lerp(t, u, v);

	q = g3[ b01 + bz1 ] ; u = ( rx0 * q[0] + ry1 * q[1] + rz1 * q[2] );
	q = g3[ b11 + bz1 ] ; v = ( rx1 * q[0] + ry1 * q[1] + rz1 * q[2] );
	b = lerp(t, u, v);

	d = lerp(sy, a, b);

	return lerp(sz, c, d);
    }

    public double noise(double vec[], int len)
    {
	switch (len)
	{
	case 0:
	    return 0.;
	case 1:
	    return noise1(vec[0]);
	case 2:
	    return noise2(vec);
	default:
	    return noise3(vec);
	}
    }

    /************************************************************
     * Methods that use the noise functions to generate height fields
     * Adapted from code written by F. Kenton Musgrave
     ************************************************************/
    /*
     * Procedural fBm evaluated at "point"; returns value stored in "value".
     *
     * Copyright 1994 F. Kenton Musgrave 
     * 
     * Parameters:
     *    ``H''  is the fractal increment parameter
     *    ``lacunarity''  is the gap between successive frequencies
     *    ``octaves''  is the number of frequencies in the fBm
     *
     * 'point' must be a double[3]
     */
    private boolean first_fBm = true;
    private double exponent_array[];
    public double
    fBm( double point[], double H, double lacunarity, double octaves )
    {
	double            value, frequency, remainder;
	int               i;
	
	/* precompute and store spectral weights */
	if ( first_fBm )
	{
	    /* seize required memory for exponent_array */
	    exponent_array = new double[(int)octaves+1];
	    frequency = 1.0;
	    for (i=0; i <= octaves; i++)
	    {
		  /* compute weight for each frequency */
		  exponent_array[i] = Math.pow( frequency, -H );
		  frequency *= lacunarity;
	    }
	    first_fBm = false;
	}
	
	value = 0.0;            /* initialize vars to proper values */
	frequency = 1.0;
	
	/* inner loop of spectral construction */
	for (i=0; i < octaves; i++)
	{
	    value += noise3( point ) * exponent_array[i];
	    point[0] *= lacunarity;
	    point[1] *= lacunarity;
	    point[2] *= lacunarity;
	}
	
	remainder = octaves - (int)octaves;
	if ( remainder != 0.0 )
	{
	    /* add in ``octaves''  remainder */
	    /* ``i''  and spatial freq. are preset in loop above */
	    value += remainder * noise3( point ) * exponent_array[i];
	}
	return( value );
    }

    /*
     * Procedural multifractal evaluated at "point"; 
     * returns value stored in "value".
     *
     * Copyright 1994 F. Kenton Musgrave 
     * 
     * Parameters:
     *    ``H''  determines the highest fractal dimension
     *    ``lacunarity''  is gap between successive frequencies
     *    ``octaves''  is the number of frequencies in the fBm
     *    ``offset''  is the zero offset, which determines multifractality
     *
     * Note: this tends to yield very small values, so the results need
     * to be scaled appropriately.
     */
    public double
    multifractal( double point[], double H, double lacunarity, 
		  double octaves, double offset )
    {
	double            value, frequency, remainder;
	int               i;

	/* precompute and store spectral weights */
	if ( first_fBm )
	{
	    /* seize required memory for exponent_array */
	    exponent_array = new double[(int)octaves+1];
	    frequency = 1.0;
	    for (i=0; i <= octaves; i++)
	    {
		  /* compute weight for each frequency */
		  exponent_array[i] = Math.pow( frequency, -H );
		  frequency *= lacunarity;
	    }
	    first_fBm = false;
	}
    
	value = 1.0;            /* initialize vars to proper values */
	frequency = 1.0;
	
	/* inner loop of multifractal construction */
	for (i=0; i < octaves; i++)
	{
	    value *= offset * frequency * noise3( point );
	    point[0] *= lacunarity;
	    point[1] *= lacunarity;
	    point[2] *= lacunarity;
	}
	
	remainder = octaves - (int)octaves;
	if ( remainder != 0.0 )
	{
	    /* add in ``octaves''  remainder */
	    /* ``i''  and spatial freq. are preset in loop above */
	    value += remainder * noise3( point ) * exponent_array[i];
	}
	return value;
    }

    /*
     * Heterogeneous procedural terrain function: stats by altitude method.
     * Evaluated at "point"; returns value stored in "value".
     *
     * Copyright 1994 F. Kenton Musgrave 
     * 
     * Parameters:
     *       ``H''  determines the fractal dimension of the roughest areas
     *       ``lacunarity''  is the gap between successive frequencies
     *       ``octaves''  is the number of frequencies in the fBm
     *       ``offset''  raises the terrain from `sea level'
     */
    public double
    Hetero_Terrain( double point[],
		double H, double lacunarity, double octaves, double offset )
    {
	double          value, increment, frequency, remainder;
	int             i;
    
	/* precompute and store spectral weights */
	if ( first_fBm )
	{
	    /* seize required memory for exponent_array */
	    exponent_array = new double[(int)octaves+1];
	    frequency = 1.0;
	    for (i=0; i <= octaves; i++)
	    {
		  /* compute weight for each frequency */
		  exponent_array[i] = Math.pow( frequency, -H );
		  frequency *= lacunarity;
	    }
	    first_fBm = false;
	}
    
	/* first unscaled octave of function; later octaves are scaled */
	value = offset + noise3( point );
	point[0] *= lacunarity;
	point[1] *= lacunarity;
	point[2] *= lacunarity;
	
	/* spectral construction inner loop, where the fractal is built */
	for (i=1; i < octaves; i++)
	{
	    /* obtain displaced noise value */
	    increment = noise3( point ) + offset;
	    /* scale amplitude appropriately for this frequency */
	    increment *= exponent_array[i];
	    /* scale increment by current `altitude' of function */
	    increment *= value;
	    /* add increment to ``value''  */
	    value += increment;
	    /* raise spatial frequency */
	    point[0] *= lacunarity;
	    point[1] *= lacunarity;
	    point[2] *= lacunarity;
	} /* for */
	
	/* take care of remainder in ``octaves''  */
	remainder = octaves - (int)octaves;
	if ( remainder != 0.0)
	{
	    /* ``i''  and spatial freq. are preset in loop above */
	    /* note that the main loop code is made shorter here */
	    /* you may want to that loop more like this */
	    increment = (noise3( point ) + offset) * exponent_array[i];
	    value += remainder * increment * value;
	}
    
	return( value );
    }
    
    
    /* Hybrid additive/multiplicative multifractal terrain model.
     *
     * Copyright 1994 F. Kenton Musgrave 
     *
     * Some good parameter values to start with:
     *
     *      H:           0.25
     *      offset:      0.7
     */
    public double 
    HybridMultifractal( double point[], double H, double lacunarity, 
			double octaves, double offset )
    {
	double          frequency, result, signal, weight, remainder; 
	int             i;
    
	/* precompute and store spectral weights */
	if ( first_fBm )
	{
	    /* seize required memory for exponent_array */
	    exponent_array = new double[(int)octaves+1];
	    frequency = 1.0;
	    for (i=0; i <= octaves; i++)
	    {
		  /* compute weight for each frequency */
		  exponent_array[i] = Math.pow( frequency, -H );
		  frequency *= lacunarity;
	    }
	    first_fBm = false;
	}
    
	/* get first octave of function */
	result = ( noise3( point ) + offset ) * exponent_array[0];
	weight = result;
	/* increase frequency */
	point[0] *= lacunarity;
	point[1] *= lacunarity;
	point[2] *= lacunarity;
	
	/* spectral construction inner loop, where the fractal is built */
	for (i=1; i < octaves; i++)
	{
	    /* prevent divergence */
	    if ( weight > 1.0 )  weight = 1.0;
	
	    /* get next higher frequency */
	    signal = ( noise3( point ) + offset ) * exponent_array[i];
	    /* add it in, weighted by previous freq's local value */
	    result += weight * signal;
	    /* update the (monotonically decreasing) weighting value */
	    /* (this is why H must specify a high fractal dimension) */
	    weight *= signal;
	
	    /* increase frequency */
	    point[0] *= lacunarity;
	    point[1] *= lacunarity;
	    point[2] *= lacunarity;
	} /* for */
	
	/* take care of remainder in ``octaves''  */
	remainder = octaves - (int)octaves;
	if ( remainder != 0.0 )
	{
	    /* ``i''  and spatial freq. are preset in loop above */
	    result += remainder * noise3( point ) * exponent_array[i];
	}
	
	return( result/2.0 - 1.0 );
    }
    
    /* Ridged multifractal terrain model.
     *
     * Copyright 1994 F. Kenton Musgrave 
     *
     * Some good parameter values to start with:
     *
     *      H:           1.0
     *      offset:      1.0
     *      gain:        2.0
     */
    public double 
    RidgedMultifractal( double point[], double H, double lacunarity,
			double octaves, double offset, double gain )
    {
	double           result, frequency, signal, weight;
	int              i;
    
	/* precompute and store spectral weights */
	if ( first_fBm )
	{
	    /* seize required memory for exponent_array */
	    exponent_array = new double[(int)octaves+1];
	    frequency = 1.0;
	    for (i=0; i <= octaves; i++)
	    {
		  /* compute weight for each frequency */
		  exponent_array[i] = Math.pow( frequency, -H );
		  frequency *= lacunarity;
	    }
	    first_fBm = false;
	}
    
	/* get first octave */
	signal = noise3( point );
	/* get absolute value of signal (this creates the ridges) */
	if ( signal < 0.0 ) signal = -signal;
	/* invert and translate (note that "offset" should be ~= 1.0) */
	signal = offset - signal;
	/* square the signal, to increase "sharpness" of ridges */
	signal *= signal;
	/* assign initial values */
	result = signal;
	weight = 1.0;
	
	for( i=1; i < octaves; i++ )
	{
	    /* increase the frequency */
	    point[0] *= lacunarity;
	    point[1] *= lacunarity;
	    point[2] *= lacunarity;
	
	    /* weight successive contributions by previous signal */
	    weight = signal * gain;
	    if ( weight > 1.0 ) weight = 1.0;
	    if ( weight < 0.0 ) weight = 0.0;
	    signal = noise3( point );
	    if ( signal < 0.0 ) signal = -signal;
	    signal = offset - signal;
	    signal *= signal;
	    /* weight the contribution */
	    signal *= weight;
	    result += signal * exponent_array[i];
	}
	
	return( (result-1.0)/2.0 );
    }

    ///** COLOR INDEX CONSTANTS ***

    public static final int BLACK = 0;
    public static final int BLUE0 = 1;
    public static final int BLUE1 = 9;
    public static final int LAND0 = 10;
    public static final int LAND1 = 18;
    public static final int WHITE = 19;

    static int[] rtable =
	{0,   0,	 0,  0,	 0,  0,	 0,  0,	 0,  0,	  0, 16, 32,
	 48, 64, 80, 96,112,128, 255};
    static int[] gtable =
	{0,   0, 16, 32, 48, 64, 80, 96,112,128, 255,240,224,208,192,
	 176,160,144,128, 255};
    static int[] btable =
	{0, 255,255,255,255,255,255,255,255,255,	  0,  4,  8,
	 12, 16, 20, 24, 28, 32, 255};

    public boolean latic; // flag for latitude based colour

    public double land;	// percentage of surface covered by land
    public double water;// percentage of surface covered by water


    ///*** METHODS REQUIRED BY INTERFACE

    /**
     * Sets internal variables required for a selected magnification,
     * image width, and image height.
     */
    public void setScaling(double M, double W, double H)
    {
    }
    /**
     * Calculates an intensity value in [0.0,1.0] at the specified point.
     */
    public double value(double x, double y, double z)
    {
	point[0] = x;
	point[1] = y;
	point[2] = z;
	switch (method)
	{
	  case METHOD_BASIC:
	    return fBm(point, H, lacunarity, octaves);
	  case METHOD_MULTIFRACTAL:
	    return multifractal(point, H, lacunarity, octaves, offset);
	  case METHOD_HETERO_TERRAIN:
	    return Hetero_Terrain(point, H, lacunarity, octaves, offset);
	  case METHOD_HYBRID_MULTIFRACTAL:
	    return HybridMultifractal(point, H, lacunarity, octaves, offset);
	  case METHOD_RIDGED_MULTIFRACTAL:
	    return RidgedMultifractal(point, H, lacunarity, octaves, offset, gain);
	}
	return 0.0;
    }
    /**
     * Returns an (alpha, red, green, blue) color value associated with
     * the value() at the specified point.
     */
    public int color(double x, double y, double z)
    {
	double alt = value(x, y, z);
	int colour;

	// calculate colour 
	if (alt <= 0.) // if below sea level then
	{
	    water++;
	    if (latic && y*y+alt >= 0.90)
	    {	// white if close to poles
		colour = WHITE;
	    }
	    else
	    {	// blue scale otherwise
		colour = BLUE1+(int)((BLUE1-BLUE0+1)*(2*alt));
		if (colour < BLUE0) colour = BLUE0;
	    }
	}
	else
	{
	    land++;
	    if (latic) alt += 0.10204*y*y;  // altitude adjusted with latitude
	    if (alt >= 0.5)	// arbitrary, but not too bad
	    {	// if high then white
		colour = WHITE;
	    }
	    else
	    {	// else green to brown scale
		colour = LAND0+(int)((LAND1-LAND0+1)*(2*alt));
	    }
	}
	if (colour < 0) colour = 0;
	if (colour > 19) colour = 19;
	return(255 << 24 | rtable[colour] << 16 | gtable[colour] << 8 | btable[colour]);
    }
    /**
     * Returns an (alpha, red, green, blue) color value associated with
     * the background value in lieu of valid noise.
     */
    public int background()
    {
	return 0xFF000000;
    }
}
